Semana 18 (Semestre 2)
So far I think I've been working blindfoled, but I did something that at least make the whole simulation finish. Before, the simulations had an error stating that there were lost atoms. And after changing so many numbers, I realized something, there was a multiplied fix ID number
And after changing this number, the simulation finished.
Al disminuar la variable Bmag de1.9428613454361806e-05 a 1.9428613454361806e-06, y dejando la simulación corrió por más tiempo, además, para esta simulación, tuve que desactivar el movimiento browniano para el Rathcet, y solo dejarlo para las particulas magnéticas.
Cabe destacar que para esta simulación, fix set_dipole fue utilizado con normalidad.
Pero se rompe al llegar al paso 7700
Pero al volver a activar brownian dynamics, la simulación únicamente llega a 7600, y de igual manera el ratchet sufre mucha deformación. Parece ser que las fluctuaciones son muy fuertes, y el ratchet debe ser reforzado de alguna forma, esta siendo agregandole una capa más, o agregando angulos entre enlaces.
Reacomodo de particulas
Para esta nueva versión, se intentará acomodar las partículas en forma de tetrahedro, para verificar si esto reduce las fluctuaciones térmicas.
Conociendo las distancias del centro del cubo a cada arista, podemos determinar los vectores base del sistema.
Podemos observar que va a ser necesario una condición que determine en que lugar del cubo estamos, por que podemos observar que en el centro de este, únicamente hay 1 partícula.
Después de modificar el código de las posiciones iniciales, por este
# Define lattice parameters
a = 0.5 # Lattice constant
radius = 10.0 # Average radius of the ratchet
saw_amp = 4.0 # Amplitude of the sawtooth modulation
saw_freq = 4 # Frequency of the sawtooth pattern (number of teeth)
N = 1000 # Number of points for the ratchet boundary
Kb = 5 #Bond energy coefficient
### For this it will be always 3D
a1 = np.array([a, 0, 0]) ##x vector
a2 = np.array([0, a, 0]) ##y vector
a3 = np.array([0, 0, a]) ##z vector
a4 = np.array([a/2, a/2, a/2])
# Number of unit cells in each direction, centered around the origin
nx, ny, nz = 3, 1, 1
# Center of the circle
center = np.array([0.0, 0.0, 0.0])
positions = []
limx, limy, limz = range(nx), range(ny), range(nz)
if nx < 2:
limx = range(nx + 1)
if ny < 2:
limy = range(ny + 1)
if nz < 2:
limz = range(nz + 1)
for i in limx:
for j in limy:
for k in limz:
pos = i * a1 + j * a2 + k * a3
positions.append(pos)
if k > 0 and j > 0 and i > 0:
positions.append([i * a/2, j * a/2, k * a/2])
positions = np.array(positions)
se obtuvo la siguiente triangulación.
Problemas con el código
El código anterior tenía problemas al posicionar las partículas centrales si el sistema crecía, así que tuve una epifania en donde conte la solución, la cuál es esta:
limx, limy, limz = range(-nx + 1, nx), range(-ny + 1, ny), range(nz)
if nx < 2:
limx = range(-nx, nx + 1)
if ny < 2:
limy = range(-ny, ny + 1)
if nz < 2:
limz = range(nz + 1)
# Loop through lattice and add particles
for i in limx:
for j in limy:
for k in limz:
pos = i * a1 + j * a2 + k * a3
positions.append(pos)
if i < max(limx):
if j < max(limy):
if k < max(limz):
pos = pos + reference
positions.append(pos)
Como podemos ver, a diferencia de una sola capa, este presenta menos fluctuaciones entre las particulas, y solo muestra rotación en todo el cuerpo.